#!/usr/bin/python
# -*- coding: utf-8 -*-
from cvxopt import matrix, spmatrix
from cvxopt.blas import gemv
from cvxopt.solvers import qp


#\begin{array}{ll} \mbox{minimize} & 2x_1^2 + x_2^2 + x_1 x_2 + x_1 + x_2 \\ \mbox{subject to} & x_1 \geq 0 \\ & x_2 \geq 0 \\ & x_1 + x_2 = 1 \end{array}
#http://abel.ee.ucla.edu/cvxopt/examples/tutorial/qp.html
Q = 2*matrix([ [2, .5], [.5, 1] ])
p = matrix([1.0, 1.0])
G = matrix([[-1.0,0.0],[0.0,-1.0]])
h = matrix([0.0,0.0])
A = matrix([1.0, 1.0], (1,2))
b = matrix(1.0)
sol=qp(Q, p, G, h, A, b)

def append_matrix_at_bottom(A, B):
    l = []
    for x in xrange(A.size[1]):
        for i in xrange(A.size[0]):
            l.append(A[i + x * A.size[0]])
        for i in xrange(B.size[0]):
            l.append(B[i + x * B.size[0]])
    return matrix(l, (A.size[0] + B.size[0], A.size[1]))


M = matrix([[4.0, 6, -4, 1.0], [6, 1, 1.0, 2.0], [-4, 1.0, 2.5, -2.0],
           [1.0, 2.0, -2.0, 1.0]])
q = matrix([12, -10, -7.0, 3])

I = spmatrix(1.0, range(M.size[0]), range(M.size[1]))
G = append_matrix_at_bottom(-M, -I)  # inequality constraint G z <= h

h = matrix([x for x in q] + [0.0 for _x in range(M.size[0])])

sol = qp(2.0 * M, q, G, h)  # find z, w, so that w = M z + q
if sol['status'] == 'optimal':
    z = sol['x']
    w = matrix(q)
    gemv(M, z, w, alpha=1.0, beta=1.0)  # w = M z + q
    print z
    print w
else:
    print 'failed'
#http://en.wikipedia.org/wiki/Linear_complementarity_problem
